I'm gonna overwrite a lot of this notebook's old content. I changed the way I'm calculating wt, and wanna test that my training worked.


In [23]:
from pearce.emulator import *
from pearce.mocks import cat_dict
import numpy as np
from os import path

In [24]:
import matplotlib
#matplotlib.use('Agg')
from matplotlib import pyplot as plt
%matplotlib inline
import seaborn as sns
sns.set()

In [25]:
from GPy.models import GPKroneckerGaussianRegression
from GPy.kern import *
import h5py

In [26]:
from pearce.emulator import GPKroneckerGaussianRegressionVar

In [27]:
training_file = '/home/users/swmclau2/scratch/xi_gg_zheng07_cosmo_v3/PearceXiggCosmo.hdf5'
test_file = '/home/users/swmclau2/scratch/xi_gg_zheng07_cosmo_test_v3/PearceXiggCosmoTest.hdf5'

In [ ]:
accs = []
log_accs = []
noise = []
for rbin in xrange(18):

    f = h5py.File(training_file, 'r')

    Ys = []
    Yerrs = []
    for i in xrange(40):
        Ys.append(f['cosmo_no_%02d'%i]['a_1.000']['obs'].value[:, rbin])
        Yerrs.append(f['cosmo_no_%02d'%i]['a_1.000']['cov'].value[:, rbin, rbin]) 

    n_hods = 100
    start_idx = 0
    X1 = f.attrs['cosmo_param_vals']
    X2 = f.attrs['hod_param_vals'][start_idx:start_idx+n_hods]
    Y = np.vstack(Ys)[:, start_idx:start_idx+n_hods]
    Yerr = np.vstack(Yerrs)[:, start_idx:start_idx+n_hods]
    
    noise.append([Y, Yerr])

    f.close()
    # how to add training errors?

    K1 =RBF(input_dim=7, ARD = False)# + RBF(input_dim=7, ARD = False)#+ Linear(input_dim = 7, ARD = False) + Bias(input_dim=7)# + White(input_dim=7)
    K2 = RBF(input_dim=4, ARD = True)# + RBF(input_dim=4, ARD = False)#+ Linear(input_dim = 4, ARD = False) + Bias(input_dim=4)# + White(input_dim=4)

    model = GPKroneckerGaussianRegressionVar(X1, X2, Y, Yerr, K1, K2)#, noise_var = 0.01 )

    model.optimize_restarts(num_restarts=10, verbose = True)

    print K1.param_array
    print K2.param_array
    print 
    
    f2 = h5py.File(test_file, 'r')

    Y2s = []
    for i in xrange(35):
        Y2s.append(f2['cosmo_no_%02d'%i]['a_1.000']['obs'].value[:, rbin])

    testX1 = f2.attrs['cosmo_param_vals']
    testX2 = f2.attrs['hod_param_vals']#[:100]
    testY = np.vstack(Y2s)#[:, :100]
    f2.close()

    predY, _ = model.predict(testX1, testX2)
    med_acc, mean_acc = np.median( np.abs( (10**predY[:,0] - 10**testY.flatten(order='F'))/10**testY.flatten(order='F') )  ), \
np.mean( np.abs( (10**predY[:,0] - 10**testY.flatten(order='F'))/10**testY.flatten(order='F') )  ) 
    print rbin, med_acc, mean_acc
    accs.append((med_acc, mean_acc))
    log_accs.append((np.median( np.abs( (predY[:,0] - testY.flatten(order='F'))/testY.flatten(order='F') )  ), \
np.mean( np.abs( (predY[:,0] - testY.flatten(order='F'))/testY.flatten(order='F') )  ) ))
    print


Optimization restart 1/10, f = -12457.688126
Optimization restart 2/10, f = -12537.0855969
Optimization restart 3/10, f = -12501.4506292
Optimization restart 4/10, f = -12566.4494755
Optimization restart 5/10, f = -12435.2716288
Optimization restart 6/10, f = -12498.8169023
Optimization restart 7/10, f = -12484.0924991
Optimization restart 8/10, f = -12494.3221698
Optimization restart 9/10, f = -12496.1188084
Optimization restart 10/10, f = -12509.1425
[ 5.6786562   3.87988861]
[ 5.76537255  1.9063896   0.68449335  2.00579988  2.55791227]

0 0.0473102786902 0.0690457923812

Optimization restart 1/10, f = -13087.9693754
Optimization restart 2/10, f = -12978.4467835
Optimization restart 3/10, f = -13079.8284725
Optimization restart 4/10, f = -13157.1642855
Optimization restart 5/10, f = -13085.8563531
Optimization restart 6/10, f = -13117.4313165
Optimization restart 7/10, f = -13083.1004805
Optimization restart 8/10, f = -13130.6212697
Optimization restart 9/10, f = -13063.5703509
Optimization restart 10/10, f = -13071.8490987
[ 15.19250025   4.89725609]
[ 13.83904354   2.26474982   0.71895658   2.50120791   6.21636794]

1 0.022749784407 0.0382069781499

Optimization restart 1/10, f = -13533.5598235
Optimization restart 2/10, f = -13661.4445731
Optimization restart 3/10, f = -13636.8328723
Optimization restart 4/10, f = -13629.880735
Optimization restart 5/10, f = -13636.4076964
Optimization restart 6/10, f = -13643.0289664
Optimization restart 7/10, f = -13650.8739466
Optimization restart 8/10, f = -13645.511221
Optimization restart 9/10, f = -13618.3162557
Optimization restart 10/10, f = -13558.2625579
[ 6.90278677  4.0065251 ]
[ 6.91287152  2.0317279   0.63224059  2.20545651  2.95543234]

2 0.0294319606301 0.0481316365829

Optimization restart 1/10, f = -13854.2476055
Optimization restart 2/10, f = -14034.4811906
Optimization restart 3/10, f = -13881.2246229
Optimization restart 4/10, f = -14073.8233192
Optimization restart 5/10, f = -14086.9594634
Optimization restart 6/10, f = -14100.6307016
Optimization restart 7/10, f = -13973.3133096
Optimization restart 8/10, f = -13979.2481077
Optimization restart 9/10, f = -13928.558213
Optimization restart 10/10, f = -13945.2793142
[ 8.74452615  4.30130321]
[ 9.88447638  2.31910481  0.67935808  2.49614309  3.33168624]

3 0.0194708492302 0.0343097936396

Optimization restart 1/10, f = -14335.3144375
Optimization restart 2/10, f = -14333.4465607
Optimization restart 3/10, f = -14331.2681037
Optimization restart 4/10, f = -14253.1504803
Optimization restart 5/10, f = -14318.964644
Optimization restart 6/10, f = -14335.2322544
Optimization restart 7/10, f = -14330.1353388
Optimization restart 8/10, f = -14325.9403671
Optimization restart 9/10, f = -14367.9970811
Optimization restart 10/10, f = -14348.0650155
[ 6.78123037  3.94502394]
[ 6.76440505  2.1395258   0.60502036  2.24698723  2.68159732]

4 0.0251307708409 0.0388964336412

Optimization restart 1/10, f = -14651.3487636
Optimization restart 2/10, f = -14664.2528455
Optimization restart 3/10, f = -14617.41534

In [ ]:
accs = np.array(accs)

In [ ]:
noise = np.array(noise)
noise.shape

In [ ]:
Y, Yvar = noise[:, 0], noise[:,1]
exp_Yvar = Yvar*np.exp(np.log(10)*Y*2) exp_errs = np.sqrt(exp_Yvar.mean(axis = (1,2)))/(np.exp(np.log(10)*Y.mean(axis = (1,2))))

In [ ]:
exp_Yerr = np.sqrt(Yvar)*10**(Y)

exp_accs = np.mean(exp_Yerr/(10**Y), axis = (1,2) )

In [ ]:
f = h5py.File(training_file, 'r')
#print f.attrs.keys()
r_bins = f.attrs['scale_bins']
f.close()
rpoints = (r_bins[1:]+r_bins[:-1])/2.0

In [ ]:
plt.plot(rpoints, exp_accs)
plt.loglog()

In [ ]:
plt.plot(rpoints, accs[:,0], label = 'Median Error')
plt.plot(rpoints, accs[:,1], label = 'Mean Error')
plt.plot(rpoints, exp_accs, label = 'Shot noise')
plt.ylim([0.001, 0.5])
plt.ylabel('Perc. Accuracy')
plt.xlabel(r'$r$ [Mpc]')
plt.legend(loc = 'best')
plt.loglog();

In [ ]:
log_accs = np.array(log_accs)

In [ ]:
plt.title('Log Acc')
plt.plot(rpoints, log_accs[:,0], label = 'Median Error')
plt.plot(rpoints, log_accs[:,1], label = 'Mean Error')
plt.plot(rpoints, (np.sqrt(Yvar)/np.abs(Y)).mean(axis = (1,2)), label = 'Shot noise')
plt.ylim([0.001, 0.5])
plt.ylabel('Perc. Accuracy')
plt.xlabel(r'$r$ [Mpc]')
plt.legend(loc = 'best')
plt.loglog();

In [ ]: